Data Analytics with R - Vendex Challenge 2 Report

Preparation

0. Import libraries and datasets

library(data.table)
library(pROC)

## Type 'citation("pROC")' for a citation.

## 
## Attaching package: 'pROC'

## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

library(stats)
library(tidyverse)

## ─ Attaching packages ──────────────────────────── tidyverse 1.3.0 ─

## ✓ ggplot2 3.2.1     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.4
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0

## ─ Conflicts ───────────────────────────── tidyverse_conflicts() ─
## x dplyr::between()   masks data.table::between()
## x dplyr::filter()    masks stats::filter()
## x dplyr::first()     masks data.table::first()
## x dplyr::lag()       masks stats::lag()
## x dplyr::last()      masks data.table::last()
## x purrr::transpose() masks data.table::transpose()

options("scipen" = 100, "digits" = 4)

machine_data = fread('/Users/robinho/Nutstore Files/01 ESADE/Study/05 MSc Business Analytics 19_20/04 Term2/01 R language/Data Analytics with R/Challenge/Challenge 2/DATA/machine_data.csv')
product_data = fread('/Users/robinho/Nutstore Files/01 ESADE/Study/05 MSc Business Analytics 19_20/04 Term2/01 R language/Data Analytics with R/Challenge/Challenge 2/DATA/product_data.csv')
transactional_data = fread('/Users/robinho/Nutstore Files/01 ESADE/Study/05 MSc Business Analytics 19_20/04 Term2/01 R language/Data Analytics with R/Challenge/Challenge 2/DATA/transactional_data.csv')
machine_failures = fread('/Users/robinho/Nutstore Files/01 ESADE/Study/05 MSc Business Analytics 19_20/04 Term2/01 R language/Data Analytics with R/Challenge/Challenge 2/DATA/machine_failures.csv')

I. Creating the relevant variable of the model

1. Merge the transactional dataset with the machine failures data set setting failure variable to 0 when no failure is recorded.

# 1-1. Merge machine_faliures with transactional data
complete_data = merge(transactional_data,
                      machine_failures,
                      by=c("machine", "timestamp", "column"),
                      all.x=T)

# 1-2. Set all the cols in "failure" for days without failure
complete_data[is.na(failure), failure:=0]
complete_data

2. In the transactional data table, create a variable called last_vend containing the timestamp of the previous sale of each machine.

# 2-1. Order the table
complete_data[order(machine, timestamp)]

# 2-2. Shift the column 'timestamp' to create the column 'last_vend'
complete_data[, last_vend := shift(timestamp,1), by = machine]
complete_data

3. Create a new variable in the transactional data table called deltahours containing, for every sale, the hours that passed since the last sale.

complete_data[, deltahours := difftime(timestamp, last_vend, units='hours')]
complete_data

4. Create an auxiliary data table called machine_daily_average with the average daily sales per machine. Use this auxiliary table to attach to every row of the transactional data table the the average daily sales per machine. You can do this by doing a merge.

# 4-1. Sales per machine
    # Comments: To calculate the transactions per machine, we consider active days to be also the ones where the machine have a failure.
num_prod = complete_data[, .(num_prod = .N), by = "machine"]
head(num_prod)

# 4-2. Number of days that each machine was working
working_days = complete_data[, .(num_days = uniqueN(date)), by = "machine"]
head(working_days)

# 4-3. Average daily sales per machine
machine_daily_average = working_days[,
                                     .(daily_sales_machine = num_prod$num_prod/num_days,
                                       machine)]

# 4-4. Get the complete_data with machine_daily_average column
complete_data = merge(complete_data,
                      machine_daily_average,
                      by = "machine")
complete_data

5. Create a new variable called delta in the transactional data table containing a normalized version of deltahours consisting on the deltahours associated with each sale divided by the average deltahours of each machine i.e. delta = deltahours /(24/daily_sales_machine). The interpretation of delta is the amount of missed sales if the machine was selling at a constant rate.

complete_data[, delta := as.numeric(deltahours) / (24 / daily_sales_machine)]

Manipulated dataset

complete_data

II. Creating the model

6. Select 30% of the machines in the transactional data for testing and 70% of the machines for training and train a linear logistic regression model called m to predict whether a machine has a failure as a function of variable delta.

# 6-1. Set the machine list
machine_list = unique(complete_data$machine)

# 6-2. Extract the training machine list & dataset
set.seed(5)
training_machine_list = sample(machine_list,
                               round(0.7*length(machine_list), 0),
                               replace = F)
train = complete_data[machine %in% training_machine_list]

# 6-3. Extract the test machine list & dataset
test_machine_list = setdiff(machine_list, training_machine_list)
test = complete_data[machine %in% test_machine_list]

# 6-4. Build the model
m = glm(failure ~ delta, data = train, family="binomial")

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

What is the value of the intercept and the coefficient accompanying variable delta?

# 6-5. Model information
summary(m)

## 
## Call:
## glm(formula = failure ~ delta, family = "binomial", data = train)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -3.037  -0.062  -0.048  -0.045   3.718  
## 
## Coefficients:
##             Estimate Std. Error z value            Pr(>|z|)    
## (Intercept) -6.91221    0.02307    -300 <0.0000000000000002 ***
## delta        0.56273    0.00314     179 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 103916  on 1287048  degrees of freedom
## Residual deviance:  58491  on 1287047  degrees of freedom
##   (1746 observations deleted due to missingness)
## AIC: 58495
## 
## Number of Fisher Scoring iterations: 9

• Comments: 1) Value of intercept: -6.91221; 2) Value of delta: 0.56273.

III. Questions related to the model

a. What’s the AUC, a measure of quality, of the model you have built on the train set? and on test set?

# a-1. Make prediction on the training dataset
pred_train0 = predict(m, newdata = train, type = "response")

# a-2. Make prediction on the test dataset
pred_test0 = predict(m, newdata = test, type = "response")

# a-3. AUC for training set
auc(train$failure, pred_train0) 

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## Area under the curve: 0.921

# a-4. AUC for test set
auc(test$failure, pred_test0)

## Setting levels: control = 0, case = 1

## Setting direction: controls < cases

## Area under the curve: 0.918

b. Plot the function of probability of failure with respect to delta to gain intuition.

# b-1. Overview of the range of "delta"
summary(complete_data$delta)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0     0.0     0.3     1.0     1.2   610.1    2495

• Comments: delta ∈ [0, 610.1]. We plot the following table, gradually diminishing the range of delta (or x), from [0, 610.1], in the end to [0, 25] because there lies the most steep part (major change of probability of failure) of the function.

# b-2. Plot the curve
curve(1/(1+exp(-(m$coefficients[[1]] + m$coefficients[[2]]*x))),
      xlim = c(0, 25),
      main = "Function of Probability of Failure",
      xlab = "Delta",
      ylab = "Probability of Failure",
      col = "#BE0606")

c. Let us create alarms with two levels of priority: med-risk and high-risk. Med-risk alarms will fire when the probability of failure is >=60% and High-risk when that probability is >=80%.

c-i. What are the threshold deltas for each type of alarm to fire?

# c-i-1. Threshold delta for medium risk alarm - "delta_0.6"
f_0.6 <- function(x) {(1/(1 + exp(-(m$coefficients[1] + m$coefficients[2]*x)))) - 0.6}
root_0.6 <- uniroot(f_0.6, c(0, 20))
delta_0.6 <- root_0.6$root
delta_0.6

## [1] 13

# c-i-2. Threshold delta for high risk alarm - "delta_0.8"
f_0.8 <- function(x) {(1/(1 + exp(-(m$coefficients[1] + m$coefficients[2]*x)))) - 0.8}
root_0.8 <- uniroot(f_0.8, c(0, 20))
delta_0.8 <- root_0.8$root
delta_0.8

## [1] 14.75

c-ii. How many of these alarms would be fired per day on average according to your model?

# c-ii-1. Calculate the total number of alarms and days (respectively) for each alarm

# Medium-risk :: total number of alarms
alarm_60 <- complete_data[delta >= delta_0.6]
nrow(alarm_60)

## [1] 3728

# Medium-risk :: total number of days
uniqueN(alarm_60$date)  

## [1] 90

# High-risk :: total number of alarms
alarm_80 <- complete_data[delta >= delta_0.8]
nrow(alarm_80)  

## [1] 2503

# High-risk :: total number of days
uniqueN(alarm_80$date)  

## [1] 90

# c-ii-2. Calculate the average number of alarms per day for each alarm

# If we calculate alarms per day for: failure_prob ∈ [0.6, +∞)
alarm_60_perDay <- nrow(alarm_60) / 90
alarm_60_perDay

## [1] 41.42

# If we calculate alarms per day for: failure_prob ∈ [0.8, +∞)
alarm_80_perDay <- nrow(alarm_80) / 90
alarm_80_perDay  

## [1] 27.81

# If we calculate alarms per day for: failure_prob ∈ [0.6, 0.8), which is not meaningful, we will not discuss this situation in future questions.
alarm_60_80_perDay <- (nrow(alarm_60) - nrow(alarm_80)) / 90
alarm_60_80_perDay  

## [1] 13.61

c-iii. What % of these will be “false alarms” i.e. failure variable is equal to 0, for each level of priority?

# failure_prob threshold == 60%
paste0('False alarm percentage for threshold 60%: ', round(complete_data[delta >= delta_0.6, .N, by = failure][failure==0]$N / nrow(alarm_60) * 100, 2), '%')

## [1] "False alarm percentage for threshold 60%: 20.57%"

# failure_prob threshold == 80%
paste0('False alarm percentage for threshold 80%: ', round(complete_data[delta >= delta_0.8, .N, by = failure][failure==0]$N / nrow(alarm_80) * 100, 2), '%')

## [1] "False alarm percentage for threshold 80%: 8.99%"

d. EWS system vs. current system

In this exercise we will estimate the profit impact of our EWS system vs the current system:
# i. If we set the EWS only with the med-risk alarms, what is the annual profit we will generate vs the current system as a % of the total profit? [For simplicity, consider the total profit to be the margin per item times the number of items in the period]
# ii. And if we set the EWS only with the high-risk alarms?
# Assumptions:
# - The time to go and fix a machine is 1.5 hours after an alarm is fired
# - The cost of an operator checking if the machine is broken is 10 euros (transport + checking). You will use this to calculate the cost of “false alarms” with our system
# - Do not consider any cost of repairing the machine since this cost needs to be done also with the current system
# - Vendex reports that their current system generates 2.2 false alarms per machine and year
# - Every item sold leaves 1.7 euros of margin on average. You will use this to calculate the impact of the “saved sales” thanks to our new EWS system

# d-i. If the med-risk alarms are set.

# 1. Select transactions that would generate an alarm.
dt_alarm_med = alarm_60  

# 2. Calculate amount of hours without sales corresponding to the threshold hours.
dt_alarm_med = dt_alarm_med[, threshold_hours := delta_0.6 * 24 / daily_sales_machine]

# 3. Add 1.5 hours to the threshold hours to take into account the repair time.
dt_alarm_med = dt_alarm_med[, threshold_hours_fixed := threshold_hours + 1.5]

# 4. Convert "threshold_hours_fixed" into "delta_fixed".
dt_alarm_med = dt_alarm_med[, delta_fixed := threshold_hours_fixed * 1 / (24 / daily_sales_machine)]

# 5. Calculate "won sales".
dt_alarm_med = dt_alarm_med[, won_sales := failure * (delta - delta_fixed)]

# 6. Calculate additional revenue.
margin_preItem = 1.7
additional_rev_med = sum(dt_alarm_med$won_sales) * margin_preItem
additional_rev_med

## [1] 69048

# 7. Calculate the cost of checking false alarms.
cost_med = sum(dt_alarm_med$failure==0) * 10
cost_med

## [1] 7670

# 8. Calculate the profit increase (in percentage) over current system
current_total_profit = margin_preItem * uniqueN(transactional_data$timestamp)
profit_increase_med = (additional_rev_med - cost_med) / current_total_profit
paste0('The med-risk alarm would generate an increase in profit of: ', round(profit_increase_med * 100, 2), '% per quarter.')

## [1] "The med-risk alarm would generate an increase in profit of: 2.29% per quarter."

=====================================================

# d-ii. If the high-risk alarms are set.

# 1. Select transactions that would generate an alarm.
dt_alarm_high = alarm_80  

# 2. Calculate amount of hours without sales corresponding to the threshold hours.
dt_alarm_high = dt_alarm_high[, threshold_hours := delta_0.8 * 24 / daily_sales_machine]

# 3. Add 1.5 hours to the threshold hours to take into account the repair time.
dt_alarm_high = dt_alarm_high[, threshold_hours_fixed := threshold_hours + 1.5]

# 4. Convert "threshold_hours_fixed" into "delta_fixed".
dt_alarm_high = dt_alarm_high[, delta_fixed := threshold_hours_fixed * 1 / (24 / daily_sales_machine)]

# 5. Calculate "won sales".
dt_alarm_high = dt_alarm_high[, won_sales := failure * (delta - delta_fixed)]

# 6. Calculate additional revenue.
margin_preItem = 1.7
additional_rev_high = sum(dt_alarm_high$won_sales) * margin_preItem
additional_rev_high 

## [1] 62437

# 7. Calculate the cost of checking false alarms.
cost_high = sum(dt_alarm_high$failure == 0) * 10
cost_high

## [1] 2250

# 8. Calculate the profit increase (in percentage) over current system
current_total_profit = margin_preItem * uniqueN(transactional_data$timestamp)
profit_increase_high = (additional_rev_high - cost_high) / current_total_profit
paste0('The high-risk alarm would generate an increase in profit of: ', round(profit_increase_high * 100, 2), '% per quarter.')

## [1] "The high-risk alarm would generate an increase in profit of: 2.25% per quarter."

=====================================================

• Comments: When using the aforementioned approach (according to the hints) to calculate the difference between delta and delta_fixed, we didn’t count the repair time for the current system; that is to say, we didn’t add 1.5 hrs to the deltahours from the last vending. So purely out of curiosity and an intention of comparison, we add the repair time to the current system to check the profit increase our system might achieve.

So again:

# d-i-2. If the med-risk alarms are set.

# 1. Select transactions that would generate an alarm.
dt_alarm_med = alarm_60  

# 2. Calculate amount of hours without sales corresponding to the threshold hours.
dt_alarm_med = dt_alarm_med[, threshold_hours := delta_0.6 * 24 / daily_sales_machine]

# 3. Add 1.5 hours to the threshold hours to take into account the repair time.
dt_alarm_med = dt_alarm_med[, threshold_hours_fixed := threshold_hours + 1.5]

# 4. Convert "threshold_hours_fixed" into "delta_fixed".
dt_alarm_med = dt_alarm_med[, delta_fixed := threshold_hours_fixed * 1 / (24 / daily_sales_machine)]

# 5. Add 1.5 hours to the hours (before broken) collected by the current system to take into account the repair time and convert it into delta-dimension info.
dt_alarm_med = dt_alarm_med[, delta_current_fixed := (as.numeric(deltahours) + 1.5)*1/(24/daily_sales_machine)]

# 6. Calculate "won sales".
dt_alarm_med = dt_alarm_med[, won_sales := failure * (delta_current_fixed - delta_fixed)]

# 7. Calculate additional revenue.
margin_preItem = 1.7
additional_rev_med = sum(dt_alarm_med$won_sales) * margin_preItem
additional_rev_med

## [1] 73165

# 8. Calculate the cost of checking false alarms.
cost_med = sum(dt_alarm_med$failure == 0) * 10
cost_med 

## [1] 7670

# 9. Calculate the profit increase (in percentage) over current system
current_total_profit = margin_preItem * uniqueN(transactional_data$timestamp)
profit_increase_med = (additional_rev_med - cost_med) / current_total_profit
paste0('The med-risk alarm would generate an increase in profit of ', round(profit_increase_med * 100, 2), '% per quarter, if we take into consideration the repair time when calculating the delta (number of sales could have taken place) of the current machine.')

## [1] "The med-risk alarm would generate an increase in profit of 2.44% per quarter, if we take into consideration the repair time when calculating the delta (number of sales could have taken place) of the current machine."

=====================================================

# d-ii-2. If the high-risk alarms are set.

# 1. Select transactions that would generate an alarm.
dt_alarm_high = alarm_80  

# 2. Calculate amount of hours without sales corresponding to the threshold hours.
dt_alarm_high = dt_alarm_high[, threshold_hours := delta_0.8 * 24 / daily_sales_machine]

# 3. Add 1.5 hours to the threshold hours to take into account the repair time.
dt_alarm_high = dt_alarm_high[, threshold_hours_fixed:= threshold_hours + 1.5]

# 4. Convert "threshold_hours_fixed" into "delta_fixed".
dt_alarm_high = dt_alarm_high[, delta_fixed := threshold_hours_fixed * 1 / (24 / daily_sales_machine)]

# 5. Add 1.5 hours to the hours (before broken) collected by the current system to take into account the repair time and convert it into delta-dimension info.
dt_alarm_high = dt_alarm_high[, delta_current_fixed := (as.numeric(deltahours) + 1.5) * 1 / (24 / daily_sales_machine)]

# 6. Calculate "won sales".
dt_alarm_high = dt_alarm_high[, won_sales := failure * (delta_current_fixed - delta_fixed)]

# 7. Calculate additional revenue.
margin_preItem = 1.7
additional_rev_high = sum(dt_alarm_high$won_sales) * margin_preItem
additional_rev_high 

## [1] 65468

# 8. Calculate the cost of checking false alarms.
cost_high = sum(dt_alarm_high$failure==0) * 10
cost_high 

## [1] 2250

# 9. Calculate the profit increase (in percentage) over current system
current_total_profit = margin_preItem * uniqueN(transactional_data$timestamp)
profit_increase_high = (additional_rev_high - cost_high) / current_total_profit
paste0('The high-risk alarm would generate an increase in profit of ', round(profit_increase_high * 100, 2), '% per quarter, if we take into consideration the repair time when calculating the delta (number of sales could have taken place) of the current machine.')

## [1] "The high-risk alarm would generate an increase in profit of 2.36% per quarter, if we take into consideration the repair time when calculating the delta (number of sales could have taken place) of the current machine."

Conclusion

• The absolute value of the increase of profit varies, depending on whether we consider repair time as a cost for the current system, which is not true if we do the calculation following the hints.
  
• In general, the med-risk alarm generates a higher increase of profit than the high-risk alarm does, which means although a threshold of 60% will lower the accuracy of predicting failure, the increase of sales (or margin) still surpasses the extra cost in transport for checking the false positive machines (FAILURE: predict-1, actual-0).